{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is based on the <a href=\"https://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html\">Comparing different clustering algorithms on toy datasets</a> example from the scikit-learn documentation.  Where possible, it has been modified to use GapStatClustering to estimate the number of custers rather than accepting the number of clusters as a given input parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(__doc__)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import time\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn import cluster, datasets, mixture\n",
    "from sklearn.neighbors import kneighbors_graph\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from itertools import cycle, islice\n",
    "\n",
    "from gapstat import GapStatClustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {\n",
    "    \"noisy_circles\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}\n",
    "    },\n",
    "    \"noisy_moons\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}\n",
    "    },\n",
    "    \"blobs\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}\n",
    "    },\n",
    "    \"no_structure\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}\n",
    "    },\n",
    "    \"aniso\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}\n",
    "    },\n",
    "    \"varied\": {\n",
    "        \"data\": None,\n",
    "        \"params\": None,\n",
    "        \"times\": {},\n",
    "        \"results\": {}        \n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============\n",
    "# Generate datasets. We choose the size big enough to see the scalability\n",
    "# of the algorithms, but not too big to avoid too long running times\n",
    "# ============\n",
    "np.random.seed(0)\n",
    "\n",
    "n_samples = 1500\n",
    "\n",
    "data[\"noisy_circles\"][\"data\"] = datasets.make_circles(n_samples=n_samples, factor=.5, noise=.05)\n",
    "data[\"noisy_circles\"][\"params\"] = {'damping': .77, 'preference': -240,\n",
    "                                   'quantile': .2, 'n_clusters': 2,\n",
    "                                   'min_samples': 20, 'xi': 0.25}\n",
    "\n",
    "data[\"noisy_moons\"][\"data\"] = datasets.make_moons(n_samples=n_samples, noise=.05)\n",
    "data[\"noisy_moons\"][\"params\"] = {'damping': .75, 'preference': -220, 'n_clusters': 2}\n",
    "\n",
    "data[\"blobs\"][\"data\"] = datasets.make_blobs(n_samples=n_samples, random_state=8)\n",
    "data[\"blobs\"][\"params\"] = {}\n",
    "\n",
    "data[\"no_structure\"][\"data\"] = np.random.rand(n_samples, 2), None\n",
    "data[\"no_structure\"][\"params\"] = {}\n",
    "\n",
    "# Anisotropicly distributed data\n",
    "random_state = 170\n",
    "X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "transformation = [[0.6, -0.6], [-0.4, 0.8]]\n",
    "X_aniso = np.dot(X, transformation)\n",
    "data[\"aniso\"][\"data\"] = (X_aniso, y)\n",
    "data[\"aniso\"][\"params\"] = {'eps': .15, 'n_neighbors': 2,\n",
    "                           'min_samples': 20, 'xi': 0.1, 'min_cluster_size': .2}\n",
    "\n",
    "# blobs with varied variances\n",
    "data[\"varied\"][\"data\"] = datasets.make_blobs(n_samples=n_samples,\n",
    "                             cluster_std=[1.0, 2.5, 0.5],\n",
    "                             random_state=random_state)\n",
    "data[\"varied\"][\"params\"] = {'eps': .18, 'n_neighbors': 2,\n",
    "                            'min_samples': 5, 'xi': 0.035, 'min_cluster_size': .2}\n",
    "\n",
    "# ============\n",
    "# Default cluster parameters\n",
    "# ============\n",
    "default_base = {'quantile': .3,\n",
    "                'eps': .3,\n",
    "                'damping': .9,\n",
    "                'preference': -200,\n",
    "                'n_neighbors': 10,\n",
    "                'n_clusters': 3,\n",
    "                'min_samples': 20,\n",
    "                'xi': 0.05,\n",
    "                'min_cluster_size': 0.1}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============\n",
    "# Cluster the datasets using different clustering algorithms\n",
    "# - use GapStatClustering when possible\n",
    "# ============\n",
    "for dataset_name in data:\n",
    "    print(dataset_name)\n",
    "\n",
    "    # update parameters with dataset-specific values\n",
    "    params = default_base.copy()\n",
    "    params.update(data[dataset_name][\"params\"])\n",
    "\n",
    "    X, y = data[dataset_name][\"data\"]\n",
    "\n",
    "    # normalize dataset for easier parameter selection\n",
    "    X = StandardScaler().fit_transform(X)\n",
    "\n",
    "    # estimate bandwidth for mean shift\n",
    "    bandwidth = cluster.estimate_bandwidth(X, quantile=params['quantile'])\n",
    "\n",
    "    # connectivity matrix for structured Ward\n",
    "    connectivity = kneighbors_graph(\n",
    "        X, n_neighbors=params['n_neighbors'], include_self=False)\n",
    "    # make connectivity symmetric\n",
    "    connectivity = 0.5 * (connectivity + connectivity.T)\n",
    "\n",
    "    # ============\n",
    "    # Create cluster objects and wrap with GapStatClustering when possible\n",
    "    # ============\n",
    "    # MeanShift is not compatible with GapStatClustering\n",
    "    ms = cluster.MeanShift(bandwidth=bandwidth, bin_seeding=True)\n",
    "    # MiniBatchKMeans without GapStatClustering\n",
    "    two_means = cluster.MiniBatchKMeans(n_clusters=params['n_clusters'])\n",
    "    # MiniBatchKMeans with GapStatClustering\n",
    "    two_means_gs = GapStatClustering(base_clusterer=cluster.MiniBatchKMeans())\n",
    "    # Ward Linkage (AgglomerativeClustering) doesn't work reliably with GapStatClustering\n",
    "    # - need to bump max_k up to a large value - need to figure out why\n",
    "    ward = cluster.AgglomerativeClustering(\n",
    "        n_clusters=params['n_clusters'], linkage='ward',\n",
    "        connectivity=connectivity)\n",
    "    # SpectralClustering without GapStatClustering \n",
    "    spectral = cluster.SpectralClustering(\n",
    "        n_clusters=params['n_clusters'], eigen_solver='arpack',\n",
    "        affinity=\"nearest_neighbors\")\n",
    "    # SpectralClustering with GapStatClustering \n",
    "    spectral_gs = GapStatClustering(base_clusterer=cluster.SpectralClustering(\n",
    "        eigen_solver='arpack',\n",
    "        affinity=\"nearest_neighbors\"))\n",
    "    # DBSCAN is not compatible wth GapStatClustering\n",
    "    dbscan = cluster.DBSCAN(eps=params['eps'])\n",
    "    # OPTICS is not compatible wth GapStatClustering\n",
    "    optics = cluster.OPTICS(min_samples=params['min_samples'],\n",
    "                            xi=params['xi'],\n",
    "                            min_cluster_size=params['min_cluster_size'])\n",
    "    # AffinityPropagation is not compatible wth GapStatClustering \n",
    "    affinity_propagation = cluster.AffinityPropagation(\n",
    "        damping=params['damping'], preference=params['preference'])\n",
    "    # Average Linkage (AgglomerativeClustering) without GapStatClustering \n",
    "    average_linkage = cluster.AgglomerativeClustering(\n",
    "        linkage=\"average\", affinity=\"cityblock\",\n",
    "        n_clusters=params['n_clusters'], connectivity=connectivity)\n",
    "    # Average Linkage (AgglomerativeClustering) with GapStatClustering \n",
    "    average_linkage_gs = GapStatClustering(base_clusterer=cluster.AgglomerativeClustering(\n",
    "        linkage=\"average\", affinity=\"cityblock\",\n",
    "        connectivity=connectivity))\n",
    "    # Birch doesn't work reliably with GapStatClustering\n",
    "    birch = cluster.Birch(n_clusters=params['n_clusters'])\n",
    "    # GaussianMixture is not currently comptible with GapStatClusterig\n",
    "    gmm = mixture.GaussianMixture(\n",
    "        n_components=params['n_clusters'], covariance_type='full')\n",
    "\n",
    "    clustering_algorithms = (\n",
    "        ('MiniBatchKMeans', two_means),\n",
    "        ('MiniBatchKMeans (GS)', two_means_gs),\n",
    "        ('AffinityPropagation', affinity_propagation),\n",
    "        ('MeanShift', ms),\n",
    "        ('Spectral', spectral),\n",
    "        ('Spectral (GS)', spectral_gs),\n",
    "        ('Ward', ward),\n",
    "        ('Agglomerative', average_linkage),\n",
    "        ('Agglomerative (GS)', average_linkage_gs),\n",
    "        ('DBSCAN', dbscan),\n",
    "        ('OPTICS', optics),\n",
    "        ('Birch', birch),\n",
    "        ('GaussianMixture', gmm)\n",
    "    )\n",
    "\n",
    "    for name, algorithm in clustering_algorithms:\n",
    "        print(\"-- \" + name)\n",
    "        t0 = time.time()\n",
    "\n",
    "        # catch warnings related to kneighbors_graph\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"the number of connected components of the \" +\n",
    "                \"connectivity matrix is [0-9]{1,2}\" +\n",
    "                \" > 1. Completing it to avoid stopping the tree early.\",\n",
    "                category=UserWarning)\n",
    "            warnings.filterwarnings(\n",
    "                \"ignore\",\n",
    "                message=\"Graph is not fully connected, spectral embedding\" +\n",
    "                \" may not work as expected.\",\n",
    "                category=UserWarning)\n",
    "            algorithm.fit(X)\n",
    "\n",
    "        t1 = time.time()\n",
    "        \n",
    "        # add run time to dataset\n",
    "        data[dataset_name][\"times\"][name] = t1 - t0\n",
    "\n",
    "        if hasattr(algorithm, 'labels_'):\n",
    "            y_pred = algorithm.labels_.astype(np.int)\n",
    "        else:\n",
    "            y_pred = algorithm.predict(X)\n",
    "        \n",
    "        # add results to dataset\n",
    "        data[dataset_name][\"results\"][name] = y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============\n",
    "# Plot the results\n",
    "# ============\n",
    "fig, axs = plt.subplots(len(data), len(clustering_algorithms), figsize=(20, 10))\n",
    "plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,\n",
    "                    hspace=.01)\n",
    "\n",
    "subplot_row = 0\n",
    "for dataset_name in data:\n",
    "    X, y = data[dataset_name][\"data\"]\n",
    "    \n",
    "    subplot_col = 0\n",
    "    for algo_name, _ in clustering_algorithms:\n",
    "        y_pred = data[dataset_name][\"results\"][algo_name]\n",
    "        runtime = data[dataset_name][\"times\"][algo_name]\n",
    "\n",
    "        if subplot_row == 0:\n",
    "            axs[subplot_row, subplot_col].set_title(algo_name, size=10)\n",
    "\n",
    "        if subplot_col == 0:\n",
    "            axs[subplot_row, subplot_col].set_ylabel(dataset_name, size=10)\n",
    "\n",
    "        colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a',\n",
    "                                             '#f781bf', '#a65628', '#984ea3',\n",
    "                                             '#999999', '#e41a1c', '#dede00']),\n",
    "                                      int(max(y_pred) + 1))))\n",
    "        # add black color for outliers (if any)\n",
    "        colors = np.append(colors, [\"#000000\"])\n",
    "\n",
    "        axs[subplot_row, subplot_col].scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])\n",
    "\n",
    "        #axs[subplot_row, subplot_col].set_xlim(-2.5, 2.5)\n",
    "        #axs[subplot_row, subplot_col].set_ylim(-2.5, 2.5)\n",
    "        axs[subplot_row, subplot_col].set_xticks(())\n",
    "        axs[subplot_row, subplot_col].set_yticks(())\n",
    "        time_text = ('%.2fs' % (runtime)).lstrip('0')\n",
    "        axs[subplot_row, subplot_col].text(0.99, 0.01, time_text, horizontalalignment='right',\n",
    "                                           size=15,\n",
    "                                           transform=axs[subplot_row, subplot_col].transAxes)\n",
    "\n",
    "        subplot_col += 1\n",
    "    subplot_row += 1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gapstat",
   "language": "python",
   "name": "gapstat"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
